/**********************************************************/
/* Aluno: Luciano G. S. Ramalho                           */
/* Numero USP: 1002101                                    */
/* Exercicio-Programa 2 -- Raizes de Equacoes Quadraticas */
/* MAC110 (BCC) -- 2007 -- Professor: Reverbel            */
/* Compilador: gcc (GCC) 4.0.3 (Ubuntu 4.0.3-1ubuntu5)    */
/**********************************************************/

/*
(a) Até que precisão funcionou o programa só com floats? 

Funcionou até epsilon = 1.0e-45.

(b) Até que precisão funcionou o programa que usa precisão
    dupla para calcular a raízes quadradas? 

Funcionou até epsilon = 1.0e-323.

(c) Até que precisão funcionou o programa que usa precisão
    extendida para calcular a raízes quadradas? 

Funcionou até epsilon = 1.0e-4950.

(d) Para cada uma das três versões do programa, explique 
    também o que acontece quando a precisão é excessiva e o 
    programa deixa de funcionar.
    
Nos três casos, quando a precisão é excessiva o laço while
da função raiz_quadrada nunca encerra, porque em nenhum 
momento o valor absoluto da diferença entre a estimativa atual 
e a anterior chega a ser menor que epsilon.

PORQUE ISSO OCORRE

Na linha 100 deste programa coloquei um assert para garantir 
que epsilon seja maior que zero, como diz o enunciado. Quando
o assert estava lá, o programa nunca travava. Rodava certo, ou 
encerrava logo na entrada com a mensagem:

  ep2-1002101.c:100: failed assertion `epsilon>0'
  Abort trap

Isso demonstra que o problema está na leitura do epsilon. 
Quando o valor fornecido para epsilon ultrapassa a precisão do 
tipo de número real escolhido, epsilon é arredondado para zero. 
Então, a expressão do laço while em raiz_quadrda sempre será 
verdadeira, pois ela só seria falsa se o valor absoluto da 
diferença fosse negativo, o que por definição não pode ocorrer.

Para entregar o exercício, comentei o assert da linha 100.

*/

#include <stdio.h>
#include <assert.h>

#define FLOAT       1
#define DOUBLE      2
#define LONG_DOUBLE 3

#define REAL DOUBLE

#if REAL == FLOAT
    typedef float real;
    #define FMT_REAL_SCAN "%f"
#elif REAL == DOUBLE
    typedef double real;
    #define FMT_REAL_SCAN "%lf"
#elif REAL == LONG_DOUBLE
    typedef long double real;
    #define FMT_REAL_SCAN "%Lf"
#endif

real absoluto(real x) {
    if (x < 0.0) return -x;
    return x;
}

real raiz_quadrada(real x, real epsilon) {
    real raiz = x;       /* estimativa */
    real raiz_ant = 0.0; /* estimativa anterior */
    while (absoluto(raiz-raiz_ant) >= epsilon) {
        raiz_ant = raiz;
        raiz = (raiz + x/raiz) / 2.0f;
    }
    return raiz;    
}

int main ()
{
    real epsilon;           /* menor valor absoluto que será considerado */
    float a, b, c;          /* coeficientes */
    float delta;            /* discriminante */
    float r1, r2;           /* raízes reais */
    float r_imag;           /* parte imaginária das raizes */
    float t;                /* varíavel temporária */
    int n;
    scanf(FMT_REAL_SCAN "%d",&epsilon,&n);
    assert(n>=1);           /* conforme o enunciado */
    /* assert(epsilon>0);*/ /* conforme o enunciado */
    printf("                    Raizes de Equacoes Quadraticas\n");
    printf("-------------------------------------------------------------------------------\n");
    printf("     coeficientes       tipo das\n");
    printf("   a      b      c      solucoes       raiz 1                raiz 2\n");
    printf("-------------------------------------------------------------------------------\n");
    for (; n>0; n--) {
        scanf("%f %f %f",&a,&b,&c);
            printf("%6.2f %6.2f %6.2f   ", a, b, c);
        if (absoluto(a) > epsilon) {
            delta = b*b-4.0*a*c;
            if (delta >= 0.0) {
                /* em vez da fórmula de Baskhara, usei uma fórmula que é 
                 * numericamente mais precisa segundo estas referências:
                 *   Introduction to Programming in Java (Sedgewick & Wayne)
                 *   http://www.cs.princeton.edu/introcs/91float/
                 *   Creative exercises #17: "Quadratic formula"
                 *   
                 *   Wikipedia
                 *   http://en.wikipedia.org/wiki/Quadratic_equation#Floating_point_implementation
                 *   Quadratic formula: Floating point implementation
                 */
                if (b > 0) {
                    t = -0.5 * (b + raiz_quadrada(delta, epsilon));
                } else {
                    t = -0.5 * (b - raiz_quadrada(delta, epsilon));
                }
                r1 = t / a;
                if (absoluto(delta) < epsilon) {
                    printf("%-14s","real dupla");
                    r2 = r1;
                } else { /* delta significativo: */
                    printf("%-14s","reais simples");
                    r2 = c / t;
                }
                /* trocar a ordem das raízes quando r1 é negativo apenas
                   para que nossa saída fique idêntica à do enunciado */
                if (r1 < 0.0) {
                    t = r1;
                    r1 = r2;
                    r2 = t;
                }
                /* "%-13s %8.4f%-13s %8.4f%s" */    
                printf("%8.4f %21.4f\n", r1, r2);
            }
            else { /* discriminante negativo */
                printf("%-14s","complexas");
                r1 = r2 = -b / (2.0*a);
                /* como as raízes complexas são conjugadas, basta calcular
                   a parte imaginária de uma delas */
                r_imag = raiz_quadrada( 4.0*a*c-b*b, epsilon) / (2.0*a);
                /* a mesma parte imaginária é usada duas vezes, com sinal trocado */
                printf("%8.4f + i*%.4f   %8.4f - i*%.4f\n", r1, r_imag, r2, r_imag);
            }
        } else {
            printf("*** ERRO: equacao nao e do segundo grau! ***\n");
        }
        
    }
    return 0;
}


